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Algols  are  the  remnants  of  binary  systems  that  have  undergone  large  scale  mass 

transfer.  This  dissertation  presents  the  results  of  the  coupling  of  a  hydrodynamical  model 

and  a  radiative  model  of  the  flow  of  gas  from  the  inner  Lagrangian  point.  The 

hydrodynamical  model  is  a  fully  Lagrangian,  three-dimensional  scheme  with  a  novel 

treatment  of  viscosity  and  an  implementation  of  the  smoothed  particle  hydrodynamics 

method  to  compute  pressure  gradients.  Viscosity  is  implemented  by  allowing  particles 

within  a  specified  interaction  length  to  share  momentum.  The  hydrodynamical  model 

includes  a  provision  for  computing  the  self-gravity  of  the  disk  material,  although  not  it  is 

not  used  in  the  present  application  to  Algols. 


VI 


Hydrogen  line  profiles  and  equivalent  widths  computed  with  a  code  by  Drake  and 
Ulrich  are  compared  with  observations  of  both  short  and  long  period  Algols.  More 
sophisticated  radiative  transfer  computations  are  done  with  the  escape  probability  code  of 
Ko  and  Kallman  which  includes  the  spectral  lines  of  thirteen  elements.  The  locations  and 
velocities  of  the  gas  particles,  and  the  viscous  heating  from  the  hydro  program  are 
supplied  to  the  radiative  transfer  program,  which  computes  the  equilibrium  temperature 
of  the  gas  and  generates  its  emission  spectrum.  Intrinsic  line  profiles  are  assumed  to  be 
delta  functions  and  are  properly  Doppler  shifted  and  summed  for  gas  particles  that  are  not 
eclipsed  by  either  star. 

Polarization  curves  are  computed  by  combining  the  hydro  program  with  the 
Wilson-Liou  polarization  program.  Although  the  results  are  preliminary,  they  show  that 
polarization  observations  show  great  promise  for  studying  circumstellar  matter. 
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CHAPTER  1 
INTRODUCTION 

The  binary  star  Algol,  also  known  as  p  Persei,  is  the  protoype  of  a  large  class  of 
eclipsing  binaries  referred  to  as  Algols.  Algols  that  eclipse  are  distinguished 
observationally  by  the  large  difference  in  depth  between  primary  and  secondary  eclipses, 
which  is  a  result  of  the  large  difference  in  temperature  between  the  two  components. 
Algols  are  also  characterized  by  emission  lines  and  light  curve  distortions  arising  from 
circumstellar  material.  The  subject  of  this  dissertation  is  the  computation  of  both  the  fluid 
flow  and  the  resulting  emission  of  radiation  by  the  gas.  By  employing  the  method  of 
smoothed-particle  hydrodynamics  (SPH),  we  compute  the  flow  of  the  gas  from  the  inner 
Lagrangian  point  (L,).  The  radiative  properties  of  the  gas  are  computed  with  the 
radiative  transfer  code  of  Ko  and  Kallman  (1994). 

In  the  early  1940's,  spectra  of  the  Algol-type  binary  RW  Tauri  taken  by  Joy  (1942) 
showed  emission  line  features  that  he  interpreted  as  arising  from  a  ring  of  circumstellar 
gas  around  the  primary  star.  Struve  observed  similar  features  in  other  Algols  (e.g.,  Struve 
1944)  and  his  immense  volume  of  high  quality  spectroscopic  observations  demonstrated 
that  these  gas  streams  and  rings,  now  more  commonly  referred  to  as  disks,  were 
characteristic  features  of  Algols  and  not  anomalies  in  just  a  few  systems.  Around  the 
same  time,  Kuiper  published  a  theoretical  treatment  of  P  Lyrae  and  pointed  out,  for  the 
first  time,  the  importance  of  the  Li  point  in  understanding  gas  flows  in  interacting 
binaries  (Kuiper,  1941). 


Spectroscopic  and  photometric  observations  of  Algols  showed  that  the  hotter  and 
more  massive  primary  stars  have  radii  that  placed  them  on  the  main  sequence  in  the  H-R 
diagram,  but  the  cooler  secondaries  were  already  evolved  from  the  main  sequence.  A  well 
known  deduction  from  stellar  evolution  theory  is  that  the  main  sequence  lifetime  of  a  star 
decreases  with  increasing  mass.  Thus,  Algols  presented  an  apparent  problem,  the  so- 
called  Algol  Paradox —  how  could  the  lower  mass  secondary  have  evolved  off  of  the 
main  sequence  before  the  higher  mass  primary  did  so?  A  major  step  in  the  resolution  of 
the  Algol  paradox  was  achieved  by  Crawford  (1955).  He  recognized  that  Algols  are  the 
remnants  of  binary  systems  that  have  experienced  large  scale  mass  transfer.  The  now  less 
massive  secondary  was  once  the  more  massive  star  and  evolved  to  fill  its  Roche  lobe, 
spilling  matter  toward  the  other  star  through  the  L,  point.  Morton  (1960)  showed  that  the 
configuration  in  which  the  more  massive  star  fills  the  Roche  lobe  is  unstable  and  results 
in  mass  transfer  on  a  thermal  time  scale.  This  meant  that  the  mass  transfer  would  be 
relatively  rapid  and  of  such  magnitude  as  to  reverse  the  mass  ratio  of  the  system,  making 
the  gainer  the  more  massive  star.  At  some  point  after  the  mass  ratio  reversal,  the  mass 
transfer  would  stabilize  and  take  place  at  the  much  longer  nuclear  time  scale.  A  key  point 
that  Morton  failed  to  take  into  account,  however,  was  the  changing  separation  of  the  stars 
as  the  mass  transfer  occured.  Later  work  by  several  researchers  allowed  for  the  changing 
separation  and  more  accurately  described  the  evolution  of  the  system  (see  the  reviews  by 
Plavec,  1968  and  Paczynski,  1971). 

Spectroscopic  observations  with  the  International  Ultraviolet  Explorer  (IUE) 
satellite  led  Plavec  (1980)  to  identify  a  class  of  binaries  whose  ultraviolet  spectra  show 
emission  lines  of  highly  ionized  species  at  all  phases,  and  he  called  these  systems  the  W 
Serpentis  stars.  The  W  Serpentis  stars  are  binaries  in  the  rapid  phase  of  mass  transfer, 
with  the  gainer  sometimes  completely  engulfed  by  the  accretion  disk,  as  is  the  case  in  P 


Lyrae  (Wilson,  1974;  Wilson  and  Terrell,  1992).  Structurally,  the  W  Serpentis  stars  are 
much  more  complicated  than  Algols  and  therefore,  much  more  difficult  to  model 
quantitatively.  One  model  of  these  systems  (Wilson,  1981)  has  the  gainer  spun  up  to  the 
centrifugal  limit  by  the  accretion  stream.  Being  at  the  centrifugal  limit,  the  gainer  cannot 
accept  any  more  of  the  incoming  gas,  and  a  thick  accretion  disk  is  formed.  Eventually  the 
gas  is  able  to  settle  onto  the  star,  and  the  system  will  then  consist  of  a  rapidly  rotating 
main  sequence  star  and  a  lobe-filling  subgiant,  making  it  a  rapidly  rotating  Algol.  If  the 
gainer  is  rotating  at  the  centrifugal  limit,  the  system  is  of  the  double-contact 
morphological  type  (Wilson,  1979).  Tidal  forces  will  cause  the  rotation  of  the  primary  to 
slow,  eventually  to  synchronism  with  the  orbit,  and  the  system  will  be  in  the  Algol  stage.1 
According  to  the  foregoing  theory,  Algols  represent  a  fairly  long-lived  state  of  certain 
binaries  that  have  undergone  large  scale  mass  transfer,  and  are  therefore  important  in 
understanding  the  evolution  of  close  binaries.  The  rapidly  rotating  Algols  are  particularly 
interesting  because  they  tend  to  be  more  active  than  the  ordinary  Algols,  but  not  quite  so 
active  as  the  W  Serpentis  systems,  which  are  more  difficult  to  model. 

The  subject  of  this  dissertation  is  the  synthesis  of  observable  quantities,  in  particular 
emission  line  profiles  and  spectral  energy  distributions,  that  arise  from  mass  transfer  in 
Algols.  We  also  present  some  results  from  the  merger  of  the  hydro  program  with  the 
Wilson-Liou  polarization  program  (Wilson  and  Liou,  1993).  Although  hydrodynamical 
treatments  and  very  simplified  radiative  treatments  of  the  problem  have  been  done 
previously,  never  before  have  they  been  done  together  so  as  to  allow  a  meaningful 
comparison  of  theory  and  observation.  The  reason  for  this  is  the  complexity  of  the 


1  This  is  true  assuming  that  the  mass-gaining  star  that  results  after  the  rapid  phase  of  mass 
transfer  is  a  main  sequence  star  which  is  what  would  be  expected  for  the  first  epoch  of 
mass  transfer.  Later  epochs  of  mass  transfer  result  in  helium  stars  or  white  dwarfs  (see 
Wilson,  1 989),  but  in  this  work  we  do  not  attempt  to  model  these  systems. 


problem  which  taxes  even  the  most  powerful  computers  available.  However,  recent 
increases  in  computing  power  combined  with  judicious  programming  and  some 
simplifying  assumptions  has  made  the  treatment  of  the  problem  feasible. 

Our  eventual  goal  is  to  have  a  computer  model  that  is  general  in  applicability  to 
matter-transferring  binaries  such  as  Algols,  cataclysmic  variables,  and  X-ray  binaries. 
Therefore,  the  implementation  of  the  model  as  a  computer  program  has  been  done  with 
generality  and  expandability  in  mind.  The  program  is  written  so  that  future  modifications 
(such  as  the  inclusion  of  magnetic  forces)  will  be  straightforward.  However,  a  generalized 
program  should  allow  the  user  to  "turn  off'  calculations  that  may  be  unnecessary  in  a 
particular  application.  An  example  occurs  in  our  present  application  to  Algols  where  the 
disk  self-gravity  is  negligible,  and  the  time-consuming  self-gravity  calculations  are 
omitted. 

The  calculation  of  the  emission  of  a  moving  gas  in  the  presence  of  the  radiation 
fields  of  two  stars  is  complicated.  The  solution  of  that  problem  will  lie  in  work  well 
beyond  the  material  presented  in  this  dissertation.  However,  by  starting  simply,  and 
gradually  increasing  the  sophistication  of  our  model,  we  believe  that  it  will  be  possible  to 
model  the  observations  of  many  types  of  interacting  binaries  accurately.  We  have  chosen 
to  begin  our  modelling  efforts  with  the  Algols  because  they  are  relatively  bright  (and 
easily  observed)  and  because  the  mass  ejection  mechanism  is  relatively  well  understood. 
Many  high-quality  observations  of  these  systems  exist.  This  dissertation  is  an  attempt  to 
bring  theory  into  comparison  with  observations. 


CHAPTER  2 
HYDRODYNAMICS 


Previous  Work 


The  earliest  attempts  to  model  the  gas  flows  in  semidetached  binaries  were 
integrations  of  the  restricted  three-body  problem  for  particles  ejected  from  the  inner 
Lagrangian  (L:)  point  at  arbitrary  angles  and  speeds.  Kopal  (1956)  computed  trajectories 
for  particles  ejected  at  various  speeds  from  L,  at  tangents  to  the  Roche  lobe  in  the  orbital 
plane.  Gould  (1957,  1959)  considered  ejections  from  L,  as  well  as  radial  ejections  from 
other  locations  along  the  equator  of  the  star,  and  considered  motion  only  in  the  orbital 
plane.  Both  authors  computed  trajectories  for  a  wide  variety  of  mass  ratios  and  ejection 
speeds  so  as  to  illuminate  the  gross  properties  of  these  flows. 

Plavec,  Sehnal,  and  Mikulas  (1964)  computed  ejections  from  L,  in  a  system  with 
mass  ratio  of  0.2151  in  an  attempt  to  model  the  circumstellar  gas  in  RW  Tauri.  Their 
results  showed  that  thermal  ejection  at  L,  was  not  sufficient  to  form  a  ring  of  material 
around  the  primary.  Only  ejections  at  higher  speeds  could  avoid  impacting  the  primary 
during  infall.  They  also  computed  trajectories  for  ejections  out  of  the  orbital  plane  and 
showed  that  such  trajectories  were  short-lived  and  could  not  encircle  the  primary.  In  a 
follow-up  paper,  Plavec  and  Kriz  (1965)  computed  trajectories  for  a  wider  range  of  mass 
ratios  (0.2,  0.3,  and  0.5)  and  for  ejections  from  various  locations  on  the  secondary's 
surface,  although  mainly  from  L,.  The  new  results  showed  that,  at  least  over  the  range  of 
mass  ratio  considered,  the  trajectories  were  surprisingly  independent  of  the  mass  ratio, 
and  that  ring-forming  trajectories  were  found  only  for  a  narrow  range  of  ejection  angles. 
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The  ring-forming  trajectories  were  also  short-lived,  indicating  that  the  disks  in  such 
systems  as  RW  Tauri  and  U  Cephei  must  be  transient  phenomena,  a  conclusion  supported 
by  observations. 

By  assuming  a  conservatively  low  density  of  10  particles  per  cm  ,  and  collision 
cross-sections  of  10"15  to  10"16  cm2,  Prendergast  (1960)  showed  that  the  mean-free  path  of 
the  particles  was  much  smaller  than  the  orbital  separation  of  the  stars  (1-10  km).  This 
result  showed  that  a  particle  trajectory  treatment  of  the  problem  was  insufficient,  and  that 
a  proper  treatment  should  be  hydrodynamical.  He  derived  approximate  solutions  to  the 
hydrodynamical  equations,  simplifying  the  problem  by  ignoring  pressure  gradient  terms 
in  Euler's  equation  and  assuming  hydrostatic  support  perpendicular  to  the  orbital  plane. 
Unfortunately,  the  computing  power  available  was  inadequate  for  a  more  general  attack 
on  the  problem,  and  Prendergast's  solutions  were  severely  limited  by  the  required 
simplifying  assumptions. 

A  more  generalized  treatment  was  given  by  Prendergast  and  Taam  (1974). 
Recognizing  that  the  equations  of  hydrodynamics  are  moments  of  the  Boltzmann 
equation,  they  simulated  solutions  of  the  Boltzmann  equation  rather  than  solving  a  set  of 
difference  equations.  The  gas  was  treated  as  a  collection  of  particles  whose  strong 
tendency  to  collide  gave  them  a  Maxwellian  distribution.  After  dividing  the  circumstellar 
region  into  a  grid  of  volumes  or  boxes,  particles  were  allowed  to  move  in  straight  lines  at 
constant  velocity  between  collisions.  After  moving  for  a  certain  interval  of  time  (the 
longest  possible  time  that  did  not  allow  any  particle  to  travel  more  than  one  box  length), 
the  particles  were  thermalized  and  a  new  Maxwellian  velocity  distribution  was 
constructed.  The  timestep  assignment  amounts  to  assigning  a  mean-free  path  for  the 
particles  (of  the  order  of  the  box  size),  so  that  viscosity  and  thermal  conductivity  are  then 
inherent  to  the  technique,  which  introduces  an  artificial  viscosity.  Unfortunately,  this 


built-in  viscosity  is  tied  to  the  computational  grid,  and  cannot  be  adjusted  independently 
of  the  grid  resolution.  The  results  of  their  calculations  for  a  system  with  parameters 
similar  to  those  of  U  Cephei  showed  that  a  stream  was  formed  that  impacted  the  primary. 
However,  some  material  did  encircle  the  primary  and  some  material  was  ejected  from  the 
system  as  a  whole,  with  more  material  being  ejected  if  the  primary  were  rotating  five 
times  faster  than  synchronously  rather  than  synchronously. 

Lubow  and  Shu  (1975)  developed  a  hydrodynamical  treatment  of  semi-detached 
binaries  that  used  matched  asymptotic  expansions.  Their  work  showed  that  particle 
trajectory  calculations  could  be  useful  over  a  large  region  that  they  called  called  the  "orbit 
region"  since  the  mean-free  path  argument  is  valid  only  if  particle  trajectories  cross  one 
another.  Their  method  makes  use  of  continuum  mechanics  and  is  semi-analytic,  meaning 
that  most  of  the  integrations  could  be  reduced  to  a  numerical  integration  of  a  set  of 
ordinary  differential  equations.  They  exploited  the  existence  of  a  parameter,  s,  which  is 
given  by  s  =  a/(Qd),  where  a  is  the  isothermal  sound  speed  (a  =  (kT/m)  ),  Q  is  the 
orbital  angular  speed  (Q  =  27r/P),  and  d  is  the  separation  of  the  two  stars.  Assuming 
synchronous  rotation  of  the  secondary  and  isothermal  flow  under  steady  conditions,  they 
reduced  the  parameter  space  of  the  equations  to  a  single  quantity —  the  mass  ratio  of  the 
binary.  They  also  presented  a  mechanism  for  the  flow  near  Lj —  a  non-isotropic  stellar 
wind  which  reaches  sonic  speeds  in  the  neighborhood  of  L,,  and  is  throttled  into  a 
narrow  stream  of  gas  making  an  angle  with  the  line  of  centers  ranging  from  19°.5  to  28°. 4 
for  the  entire  range  of  possible  mass  ratios  (cf.  their  Table  1).  In  a  later  paper  (Lubow  and 
Shu,  1 976),  they  considered  the  dynamics  of  the  gas  flow  perpendicular  to  the  orbital 
plane,  and  found  that  the  scale  height  of  the  stream  often  exceeded  the  corresponding 
hydrostatic  value  because  the  inertia  of  the  gas  prevented  it  from  adjusting  to  the  rapidly 
varying  gravitational  field. 


About  the  same  time,  Lin  and  Pringle  (1976)  outlined  a  2-dimensional,  many-body 
approach  to  the  situation  where  the  size  of  the  primary  star  was  very  small  compared  to 
the  orbital  separation,  such  as  in  cataclysmic  variables  and  X-ray  binaries.  Their  method 
was  the  first  fully  Lagrangian  method  for  simulating  gas  flows  in  interacting  binaries,  and 
included  a  viscous  interaction  for  particles,  but  did  not  include  pressure  gradients.  The 
viscosity  was  effected  by  allocating  the  particles  to  the  boxes  of  a  Cartesian  grid  centered 
on  the  primary  star,  and  then  allowing  the  particles  within  a  box  to  interact  with  one 
another  so  that  the  particles  achieved  solid-body  rotation  about  their  center  of  mass. 
Unfortunately,  this  scheme  has  the  undesirable  consequence  of  allowing  particles  that  are 
close  to  one  another,  but  not  in  the  same  box,  to  avoid  interaction  while  more  distant 
particles  in  the  same  box  will  interact.  Lin  and  Pringle  concluded  that  the  accretion  disk 
was  quite  well  defined  and  comparable  in  size  to  the  Roche  lobe,  with  all  but  a  few 
percent  of  the  transferred  mass  being  accreted  by  the  primary. 

Whitehurst  (1988b)  presented  an  extension  of  the  Lin  and  Pringle  scheme  that 
included  pressure  gradients  computed  according  to  the  scheme  of  Larson  (1978),  which  is 
a  simplistic  approach  necessary  to  keep  computational  requirements  reasonable,  but  it  is 
not  as  drastic  as  might  first  appear  since  the  disk  is  dominated  by  angular  momentum 
transport.  Larson's  method  assumes  that  the  particles  are  extended,  deformable  gas  clouds 
in  continual  contact.  The  repulsive  acceleration  of  the  particles  is  given  as  C  /r,  where  C 
is  a  parameter  proportional  to  the  local  sound  speed  and  r  is  the  separation  of  the 
particles.  Whitehurst's  major  simplifying  assumption  was  that  the  energy  dissipated  by 
particle  interactions  was  radiated  away  instantaneously  through  the  disk  surfaces.  This,  of 
course,  made  the  pressure  calculations  somewhat  crude,  but  it  was  the  next  obvious  step 


in  simulating  flows.  Whitehurst  (1988a)  applied  his  model  to  the  SU  UMa  star  Z 
Chamaeleontis  and  was  able  to  model  the  superhumps1  of  the  system  successfully. 

The  Roche  Model 

The  basic  framework  for  our  modelling  of  Algols  is  the  extended  Roche  model, 
which  is  based  on  several  assumptions  about  the  stars.  These  assumptions  for  the 
standard  Roche  model  are  as  follows: 

1.  Gravitationally,  the  stars  act  as  point  masses.  This  assumption  is  justified 
because  of  the  high  central  condensation  of  stars.  In  Algols,  the  more 
distorted  secondary  has  evolved  beyond  the  main  sequence  and  is  even  more 
centrally  condensed. 

2.  The  stars  orbits  are  circular.  Light  curves  of  Algols  show  that  this  assumption 
is  valid  since  the  eclipses  are  equally  spaced  and  of  equal  duration. 

3.  The  stars  rotate  synchronously.  For  the  distorted  secondary,  synchronous 
rotation  is  invariably  observed  because  of  the  large  tidal  torques  produced  on 
it.  However,  the  primary  components  of  many  Algols  do  not  rotate 
synchronously,  and  the  Roche  model  must  be  extended  to  treat  them 
properly. 

With  the  above  assumptions,  one  can  use  a  frame  of  reference  that  rotates  with  the 
binary,  and  thus  has  the  stars  at  fixed  locations.  The  surfaces  of  the  stars  will  coincide 
with  surfaces  of  constant  potential  energy,  or  equipotentials.  There  exists  one 
equipotential  of  critical  importance  for  close  binaries.  That  is  the  one  that  includes  the  L, 


'Superhumps  are  periodic  increases  in  brightness  of  amplitude  20-30  percent  that  occur 
during  superoutbursts  (outbursts  lasting  10  to  14  days  as  opposed  to  normal  outbursts 
which  last  2  to  3  days). 
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point.  This  potential  defines  the  Roche  lobe  of  each  star,  beyond  which  it  cannot  expand. 
Figure  2-1  shows  the  Roche  lobes  for  the  SX  Cas  system.  At  the  L,  point,  which  lies 
along  the  line  of  centers  closer  to  the  less  massive  secondary  component,  the  sum  of  the 
gravitational  and  centrifugal  forces  vanishes.  If  one  of  the  stars  tries  to  expand  beyond  its 
Roche  lobe,  the  material  at  the  Lj  point,  with  gas  pressure  on  one  side  and  essentially  a 
vacuum  on  the  other,  will  flow  toward  the  other  star.  In  normal  Algols,  the  secondary  fills 
its  Roche  lobe  and  dumps  material  toward  the  detached  primary  intermittently  while 
expanding  on  a  nuclear  timescale. 
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Figure  2-1.  The  Roche  lobes  for  the  components  of  SX  Cas.  The  surface  of 
the  detached  primary  star  is  indicated  by  the  circle. 
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The  extended  Roche  model  allows  for  a  more  general  treatment  of  close  binaries.  In 
this  model,  the  stars  may  have  eccentric  orbits  and  rotate  non-synchronously  (see  Wilson, 
1979).  Although  eccentric  orbit  generalization  is  not  necessary  for  Algols,  non- 
synchronous  rotation  has  proven  to  be  reasonably  common  for  the  primary  components. 
In  the  rapidly  rotating  Algols,  the  primaries  do  not  have  a  Roche  lobe  per  se,  but  they  do 
have  a  critical  surface.  The  critical  surface  is  defined  by  the  equipotential  that,  for  a  given 
rotation  rate,  contains  a  point  at  the  centrifugal  limit.  Double-contact  binaries  have  the 
secondary  fdling  its  Roche  lobe  while  the  primary  fills  its  rotational  critical  surface 
(Wilson,  1979). 

Our  Scheme 

A  large  number  of  techniques  have  been  developed  to  solve  the  equations  of 

hydrodynamics,   but  they  can  be  classified  into  three  groups:   Lagrangian  (particle 

schemes),  Eulerian  (grid  or  mesh  schemes),  and  a  combination  of  the  two.  Each  technique 

has  advantages  and  disadvantages  that  must  be  carefully  considered  for  applications  to 

different  situations.  Non-adaptive  grid  codes  are  straightforward  to  program  and  analyze, 

but  can  be  very  inefficient  when  the  density  of  the  fluid  changes  rapidly.  They  require 

large  amounts  of  memory  when  applied  to  three-dimensional  flows.  Particle  methods,  on 

the  other  hand,  benefit  from  the  fact  that  the  number  density  of  the  particles  maps  directly 

onto  the  mass  density  of  the  fluid,  giving  increased  resolution  in  high  density  areas. 

Other  advantages  of  a  particle  method  make  it  attractive  for  computing  gas  flows  in 
binary  star  systems.  One  attractive  feature  is  that  there  is  no  grid  to  impose  an  artificial 
geometry  on  the  system.  Another  is  the  fact  that  a  particle  scheme  can  be  "low  level," 
which  means  that  particle  motions  are  computed  by  integrating  the  forces  acting  on  them, 
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which  makes  the  addition  of  new  effects  (e.g.,  disk  self-gravity  or  magnetic  fields) 
possible  with  minimal  coding  effort. 

Therefore,  we  have  adopted  a  three-dimensional  Lagrangian  scheme  similar  to  the 
two-dimensional  Lin  and  Pringle  (1976)  and  Whitehurst  (1988b)  approaches.  An  early 
version  of  our  program  was  described  previously  (Terrell  and  Wilson,  1993). 
Assumptions  about  the  binary  system  appropriate  to  Algols  are  made,  namely  that  the 
stars  act  gravitationally  as  centrally  condensed  mass  points  and  move  in  circular  orbits. 
Particle  positions  and  velocities  are  specified  in  a  rectangular  coordinate  system  with  the 
origin  at  the  center  of  mass  of  the  binary,  and  the  unit  of  distance  being  the  stellar 
separation.  The  x-axis  of  the  system  lies  along  the  line  of  centers  and  is  positive  in  the 
direction  of  the  lobe-filling  component.  The  z-axis  is  parallel  to  the  orbital  angular 
momentum  vector  of  the  binary,  and  the  y-axis  is  defined  such  that  the  coordinate  system 
is  right-handed.  The  unit  of  time  is  P/2%,  where  P  is  the  orbital  period.  The  mass  ratio  is 
defined  in  the  sense  of  loser  to  gainer, 


We  define  an  auxiliary  quantity 


*  =  t  (2-D 


(2-2) 


l  +  q 


and  write  the  acceleration  components  due  to  gravity  of  a  particle  located  at  (x,y,z)  as 

1  -  u  u 

x-2y  +  x- — ,—  (x-x. )-— r(x-x2)  (2-3) 

y^-lx  +  y-^y-^y  (2-4) 

z  =  -^-z-^z  (2-5) 

r\  r2 

where 

rf  =  (x-x,)2  +  v2  +  z2  (2-6) 
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with  x\  and  x2  being  the  distances  from  the  center  of  mass  of  the  primary  and  secondary, 
respectively. 

Particles  are  ejected  from  L,with  initial  velocities  specified  in  the  program  input  file. 
The  user  can  specify  an  angle  (J)  that  is  the  projection  of  the  initial  velocity  in  the  orbital 
plane,  measured  from  the  positive  x-axis.  (Thus  4>=1 80°  means  the  particle  is  ejected 
toward  the  center  of  the  primary  star.)  The  angle  8  is  measured  from  the  positive  z-axis  to 
the  ejection  vector,  so  0=90°  means  that  the  ejection  is  in  the  orbital  plane  and  0=0° 
means  that  ejection  is  parallel  to  the  orbital  angular  momentum.  Angles  A<|>  and  A6  are 
user-specified  ranges  over  which  the  particles  will  be  randomly  ejected,  centered  on  <j)  and 
0.  That  is,  the  particles  will  be  ejected  randomly  with  initial  angles  <|>  ±  A<j>/2  and  6  ± 
A9/2.  The  initial  velocity  specification  is  completed  by  the  initial  dimensionless  speed.  In 
our  experiments  the  initial  ejection  velocities  were  consistent  with  the  calculations  of 
Lubow  and  Shu  (1975). 

Of  critical  importance  in  simulating  astrophysical  disks  is  the  handling  of  viscosity. 
The  physical  mechanism  of  the  viscosity  in  accretion  disks  is  not  known,  although 
turbulence  and  magnetic  fields  have  been  suggested  (see  Shakura  and  Sunyaev  (1973) 
and  Pringle  (1981)).  We  have  chosen  to  implement  the  viscosity  by  specifying  a  viscous 
interaction  length,  a,  for  the  particles.  Particles  within  a  distance  a  of  one  another  are 
allowed  to  share  momentum  if  they  are  approaching  one  another.  At  each  timestep  all 
particles  within  a  of  a  particular  particle  are  found.  Label  the  particle  whose  new  velocity 
we  wish  to  compute  by  the  number  1,  and  the  interacting  particles  by  2,3,4...  n+\  (where 
n  is  the  number  of  particles  interacting  with  particle  1),  then  compute  the  following 
quantities: 

(xl+x2Xl-r,2)  +  (x]+i,)(l-ri3)+...+(x,+i„Xl-^) 

*„  =  -  -a-  -*-  -^      (2-8) 

n 
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x,r12  +r,3+...+  r,"+i 
^=-^— - 

n 


The  new  x-component  of  the  velocity  for  particle  1  is  then  computed  by 

xnew  =  xavg+xkeep  (2-10) 

and  the  y  and  z-components  are  computed  in  the  same  manner.  The  particles  retain  their 
positions  while  having  their  velocities  modified  in  a  way  that  conserves  momentum  and 
does  not  depend  on  the  order  of  calculation  as  in  the  method  used  in  Terrell  and  Wilson 
(1993). 

Inefficiency  in  the  determination  of  nearest  neighbors  can  cause  rapid  escalation  of 
computing  time.  The  simplest,  and  most  inefficient  method  would  be  to  compare  each  of 
n  particles  with  the  others,  requiring  54  «(«-l  Comparisons.  To  improve  the  efficiency  of 
this  process,  we  have  implemented  a  grid  location  scheme  to  eliminate  particles  too 
distant  from  a  given  particle  to  interact  with  it.  A  grid  centered  on  the  primary  star  with 
box  sizes  of  length  a  is  used  to  generate  a  linked  list  of  interacting  particles.  Starting  with 
the  first  particle,  its  grid  location  is  determined  and  the  particle  number  is  placed  in  a  2- 
dimensional  array,  GRID(x,y),  corresponding  to  the  grid.  Then  then  grid  location  of 
particle  two  is  determined.  If  particle  2  happens  to  be  in  the  same  box  as  particle  1,  then 
an  auxiliary  1 -dimensional  array  POINT  is  updated  so  that  POINT(2)  contains  a  1, 
indicating  that  particle  2  "points"  to  particle  1 .  Now  suppose  particle  3  is  also  in  the  same 
box.  The  value  in  GRID  is  changed  to  a  3,  and  POINT(3)  is  set  to  2.  This  scheme  has  the 
advantage  that  it  uses  only  n  memory  locations.  A  simpler,  but  much  less  efficient 
scheme  would  be  to  make  GRID  a  3 -dimensional  array,  with  dimensions  1  and  2  being 
the  x  and  y  locations  of  the  grid  box  and  dimension  3  being  the  labels  of  the  particles  in 
that  grid  box.  Although  very  easy  to  implement,  this  latter  method  is  very  inefficient 
because  each  grid  box  must  be  dimensioned  to  n  in  the  third  dimension  for  every  grid  box 
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in  order  to  avoid  possible  overflow  errors.  And  if  the  search  grid  must  be  extended  to 
three  dimensions  (for  geometrically  thick  disks),  GRID  would  have  to  be  a  4-dimensional 
array,  further  increasing  the  memory  requirements. 

Once  all  of  the  particles  have  been  assigned  to  the  search  grid,  interacting  particles 
can  be  quickly  identified.  To  find  all  particles  interacting  with  particle  /',  one  need  only 
search  the  box  containing  particle  i  and  the  eight  boxes  that  surround  it  (assuming  the  box 
containing  particle  i  is  not  on  the  edge  of  the  grid). 

The  viscous  interactions  between  the  particles  cause  their  total  kinetic  energy  to 
decrease,  and  this  loss  of  kinetic  energy  shows  up  in  the  form  of  thermal  energy.  The 
hydro  program  computes  the  total  viscous  heating  at  each  timestep  by  comparing  the  new 
kinetic  energy  of  each  particle  with  its  previous  value.  The  acquired  thermal  energy  is 
then  distributed  to  each  particle  based  on  the  strength  of  its  viscous  interaction.  The 
radiative  transfer  program  uses  the  viscous  heating  rate  (in  ergs/cm  /sec)  to  compute  the 
equilibrium  temperature  of  the  gas.  The  heating  rate  is  computed  by  assuming  that  the 
viscously  deposited  energy  is  added  to  the  gas  evenly  over  the  timestep,  and  evenly  over 
the  volume  of  the  gas  particle. 

Forces  arising  from  pressure  gradients  in  the  disk  are  typically  small  compared  to 
gravitational  and  centrifugal  forces,  but  can  be  important  nonetheless  since  the  latter  two 
forces  can  be  anti-parallel,  thus  may  largely  cancel  one  another  and  make  the  pressure 
gradient  forces  non-negligible.  Our  code  uses  the  smoothed  particle  hydrodynamics 
(SPH)  scheme  developed  by  Lucy  (1977)  and  Gingold  and  Monaghan  (1977,  1978, 
1 979a,  and  1 979b)  to  avoid  certain  limitations  of  grid-based  finite  difference  schemes. 
Since  its  introduction,  SPH  has  become  a  standard  tool  in  computational  astrophysics, 
with  applications  to  diverse  areas  such  as  the  formation  of  planetoids  (Benz,  et  al.,  1989), 
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star  formation  (Prongracic,  et  al.,  1993),  galactic  collisions  (Hernquist  and  Katz,  1989), 
and  the  formation  of  galactic  clusters  (Thomas  and  Couchman,  1992). 

Most  of  the  advantages  of  SPH  arise  from  the  fact  that  it  treats  fluid  elements  as 
extended  clouds  of  material  whose  centers  of  mass  move  according  to  the  conservation 
laws  of  hydrodynamics.  Unlike  grid-based  methods,  there  are  no  artificial  constraints  on 
the  geometry  of  the  system  or  the  volume  that  the  fluid  occupies.  SPH  also  has  the 
advantage  that  physical  laws  are  applied  in  their  most  fundamental  form  since  the 
interaction  between  particles  is  specified. 

To  compute  values  of  continuum  variables  such  as  pressure  and  density  from  the 
discrete  distribution  of  particles,  an  average  over  the  particle  distribution  is  taken  by 
means  of  an  appropriate  smoothing  function.  This  procedure  has  the  advantage  that 
spatial  derivatives  of  the  continuous  quantities  can  be  replaced  by  the  analytically 
calculated  derivatives  of  the  smoothing  function.  The  smoothing  function,  w(r-r'), 

must  satisfy 

\w{r-r;)d'r  =  \  (2-11) 

v 

where  V  represents  the  volume  of  interest  and  r  is  the  position  vector  of  particle  j.  With 
n  particles  of  mass  m,  the  smoothed  density  is  given  by 

n 

p(F)  =  mIvc(/:-r)  (2-12) 

and  the  acceleration  acting  on  particle  i  due  to  pressure  gradients  is 


dv 


rp jrs 


—L  =  -mY  -^  +  ^    VMfi-rj)  (2-13) 

dt  i\Pi      PjJ 

Schussler  and  Schmitt  (1981)  discuss  the  pitfalls  of  using  certain  smoothing 

functions,  and  recommend  the  following  smoothing  function  for  modelling  fluids  where 

self-gravity  is  negligible: 


w(r)  =  ' 


,  r 
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for  r  <■  h 

(2-14) 


0  for  r  >  /? 

where  the  normalization  constant  is  given  by 


21 
47t/?3 


c  =  ^  (2-15) 


and  h  is  the  user-specified  smoothing  length.  This  smoothing  function  avoids  artificial 
clustering  that  occurs  in  situations  where  self-gravity  is  negligible  (such  as  Algol  disks), 
and  we  employ  it  in  our  program. 

In  some  disk  models  such  as  for  W  Serpentis  stars,  disk  self-gravity  must  be 
considered  (See  Wilson  (1981)  and  Wilson  (1974)  for  a  discussion  of  the  disk  of  p 
Lyrae.).  Because  our  program  was  designed  to  be  extended  easily,  we  were  able  to  add 
the  computation  of  forces  due  to  disk  self-gravity  with  very  little  effort,  although  we  do 
not  use  it  in  our  present  application  to  Algols.  The  user  can  enable  the  self-gravity 
calculations  by  setting  a  control  integer  in  the  input  file.  To  reduce  the  number  of 
computations,  the  gravitational  interaction  between  particles  is  computed  with  a  hybrid 
grid  and  n-body  scheme.  Particles  close  to  one  another  are  treated  as  individual  sources  of 
gravity  and  the  forces  are  computed  in  the  usual  softened  n-body  manner.  More  distant 
particles  are  treated  as  groups  and  the  force  is  computed  based  on  the  sum  of  the  masses 
and  their  common  center  of  mass. 

Some  Algols  have  primary  components  that  rotate  faster  than  synchronously 
because  of  spin-up  during  mass  transfer.  To  simulate  these  systems,  it  is  necessary  to  to 
model  the  interaction  of  the  impacting  stream  and  the  surface  of  the  star.  In  our  model 
this  is  done  by  allowing  stream  particles  to  interact  with  surface  particles  moving  at 
velocities  specified  by  an  input  parameter  F,  which  is  defined  as  the  angular  rotation  rate 
of  the  star  divided  by  the  orbital  angular  speed  (see  Wilson,  1979).  Thus  F=\  for 
synchronous  rotation.  The  program  also  allows  for  a  "bounce"  parameter  for  specifying 
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the  strength  of  the  stream-surface  interaction.  The  radial  component  of  the  stream 
velocity  at  the  star  surface  is  reversed  and  multiplied  by  the  bounce  parameter.  With  the 
bounce  parameter  set  to  0,  the  stream  plows  into  the  star  and  is  engulfed  by  it.  A  value  of 
1  for  the  bounce  parameter  causes  a  perfectly  elastic  collision. 

Finally,  the  mass  transfer  rate  (in  M0  yr"  )  is  specified  in  the  input  file  along  with 
the  number  of  timesteps  between  the  ejection  of  particles.  Once  these  two  parameters 
have  been  chosen,  the  mass  of  each  gas  particle  is  determined.  Particles  will  be  emitted 
from  L,  until  the  user-specified  limit  has  been  reached  and  the  program  will  run  until  all 
particles  have  either  been  accreted  by  either  star  or  lost  from  the  binary  system.  (Note  that 
the  length  of  the  mass  transfer  event  is  determined  by  the  maximum  number  of  particles 
and  the  number  of  timesteps  between  ejections.) 

Numerical  Models 

Algol  systems  can  be  split  into  two  groups  in  terms  of  the  type  of  disk  formed 

during  mass  transfer.  The  difference  between  the  two  is  determined  by  whether  or  not  the 
matter  stream  strikes  the  primary  star  before  encircling  it.  Lubow  and  Shu  (1975)  in  their 
Table  2  list  a  parameter,  tnmjn,  which  is  the  minimum  distance  of  the  stream  from  the 
center  of  the  primary  star  (in  units  of  the  orbital  separation).  Figure  2-1  shows  a  plot  of 
relative  primary  radius  versus  mass  ratio  (r-q  diagram)  with  tumin  indicated  by  the  smooth 
curve  and  the  values  for  the  primaries  of  a  number  of  Algols  and  W  Serpentis  stars.  If  the 
radius  of  the  primary  star  is  larger  than  Tumin,  then  the  stream  will  impact  the  star,  and  a 
transient  disk  will  form.  If  the  radius  of  the  star  is  less  than  t3min,  then  the  stream  flows 
around  the  star  and  impacts  itself,  as  seen  in  a  simulation  of  SX  Cassiopeiae  shown  in 
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Figure  2-2.  Relative  radius  versus  mass  ratio  for  several  Algols  and  W  Serpentis 
stars.  The  smooth  curve  indicates  the  Lubow-Shu  stream  impact 
condition.  Systems  observed  to  have  permanent  disks  are  indicated  by 
triangles  and  those  having  transient  disks  by  circles. 


Figure  2-3.  There  are  some  systems  in  Figure  2-2  that  seem  not  to  obey  this  logic,  in 
particular  V356  Sgr,  S  Cnc,  and  RW  Per.  In  the  case  of  V356  Sgr,  the  primary  is  most 
likely  rotating  at  the  centrifugal  limit  (Wilson  and  Caldwell,  1978),  thus  inhibiting  the 
accretion  of  the  circumstellar  material.  Although  its  location  in  the  r-q  diagram  indicates 
that  it  should  have  a  permanent  disk,  S  Cnc  shows  very  little  activity,  which  is  not 
suprising  given  its  rather  extreme  mass  ratio.  It  does,  however,  exhibit  rapid  rotation  (Van 
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Hamme  and  Wilson,  1993).  RW  Per  was  long  thought  to  have  a  thick  accretion  disk 
(Hall,  1969;  Hall  and  Stuhlinger,  1978),  but  Wilson  and  Plavec  (1988)  found  no  evidence 
of  such  a  disk  in  their  light  curve  modelling.  They  conclude  that  RW  Per  is  an  "ordinary 
Algol  system"  with  the  rapidly  rotating  primary  filling  about  94%  of  its  limiting  lobe.  So, 
although  there  are  some  anomalies  in  the  r-q  diagram,  the  Lubow-Shu  criterion  does  seem 
to  be  a  generally  good  indicator  of  the  type  of  disk  that  will  form  in  a  binary. 

For  a  given  absolute  radius,  the  relative  radius  of  the  primary  will  scale  as  P~  , 
where  P  is  the  orbital  period.  Thus  we  expect  long-period  systems  generally  to  have 
permanent  disks  and  short-period  systems  to  have  transient  disks.  Peters  (1989)  finds  that 
systems  with  P  <  4.5  days  rarely  show  permanent  Ha  emission,  while  those  with  P  >  6 
days  generally  exhibit  permanent  Ha  emission.  Highly  variable  emission  is  exhibited  by 
systems  whose  periods  lie  in  the  range  of  5-6  days.  Computed  Ha  profiles  for  both 
transient  and  permanent  disk  systems  will  be  presented  in  Chapter  4. 

Figure  2-4  and  Figure  2-5  show  the  formation  of  the  SX  Cas  disk  as  a  mass  transfer 
event  continues.  In  these  simulations,  the  viscous  interaction  parameter  was  0.06.  By  1.6 
orbital  periods  an  elliptical  disk  has  formed.  For  comparison,  a  simulation  with 
a  =  0.1  was  also  run.  Figure  2-6  shows  the  SX  Cas  disk  at  the  same  time  as  shown  in 
Figure  2-5.  With  stronger  viscous  interaction,  the  disk  shows  greater  radial  spreading,  as 
expected,  and  is  also  elliptical.  The  elliptical  shape  is  maintained  to  at  least  8  orbits,  when 
the  simulation  was  stopped.  Figure  2-7  shows  the  disk  after  8  orbits  and  the  gravitational 
influence  of  the  secondary  is  seen  in  spiral  arms  that  have  developed.  The  spiral  arms  are 
caused  by  the  3:1  resonance  which  leads  to  stresses  in  the  disk.  This  phenomenon  is  seen 
in  the  simulations  of  Whitehurst  (1988a,  1988b)  and  Lubow  (1992)  gives  a  detailed 
discussion  of  the  instability  mechanisms. 
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Figure  2-3.  A  simulation  of  SX  Cas,  showing  the  stream  trajectory  0.38  orbital  periods 
after  the  beginning  of  the  transfer  event. 
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Figure  2-4.  The  early  stages  of  disk  formation  in  SX  Cas  at  0.8  orbital  periods  after  the 
onset  of  mass  transfer. 
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Figure  2-5.  The  elliptical  disk  that  forms  1.6  orbital  periods  after  the  onset  of  mass 
transfer  in  SX  Cas. 
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Figure  2-6.  A  simulation  at  the  same  phase  and  with  the  same  parameters  used  for  Figure 
2-5,  except  that  the  viscous  interaction  length  is  0.10  as  opposed  to  0.06. 
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Figure  2-7.  The  same  simulation  as  in  Figure  2-6  at  8  orbits  after  the  onset  of  mass 
transfer.  Note  the  spiral  arms  that  have  formed  due  to  the  3:1  resonance. 


CHAPTER  3 
RADIATIVE  TRANSFER 


The  Escape  Probability  Technique 

An  early  version  of  our  program  (Terrell  and  Wilson  ,1992)  used  the  hydrogen-only 
radiative  transfer  code  of  Drake  and  Ulrich  (1980;  hereafter  DU)  to  compute  Ha  profiles. 
This  code  is  useful  for  computing  line  profiles  because  its  computational  demands  are 
fairly  modest.  Many  ground-based  observations  of  Algols  are  available  (mainly  Ha  and 
Hp)  and  the  DU  code  provides  a  means  of  rapidly  computing  hydrogen  line  profiles.  DU 
include  all  relevant  collisional  processes,  and  explicitly  treat  the  first  six  angular 
momentum  sublevels  for  the  first  six  energy  levels.  The  angular  momentum  sublevels  for 
higher  energy  levels  are  computed  by  assuming  statistical  equilibrium.  Justifications  for 
these  assumptions  are  given  by  DU,  and  the  assumptions  are  valid  for  the  conditions 
found  in  Algols. 

Unfortunately,  the  only  hydrogen  line  observable  by  the  International  Ultraviolet 
Explorer  (IUE)  is  Lyman  a,  although  there  are  many  lines  of  other  species  in  the 
wavelength  range  covered  by  the  IUE  short  wavelength  camera  (see  Figure  3-1). 
Therefore,  we  require  the  capability  to  include  other  atomic  species  such  as  helium, 
carbon,  and  nitrogen.  Fortunately  the  necessary  program  was  completed  recently  in  the 
doctoral  dissertation  of  Yuan-Kuen  Ko  at  the  University  of  Maryland  (Ko  and  Kallman, 

24 


25 

1994).  This  program  includes  the  lines  of  the  thirteen  cosmically  abundant  elements —  H, 
He,  C,  N,  O,  Ne,  Mg,  Si,  S,  Ca,  Fe,  Ni  and  Ar.  The  following  discussion  covers  the 
general  techniques  employed  in  the  program.  For  details  the  reader  is  referred  to  Ko  and 
Kallman(1994). 
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Figure  3-1.  An  IUE  spectrum  of  the  Algol-type  binary  U  Cephei  during  its  total  primary 
eclipse  showing  the  rich  set  of  emission  lines  arising  from  the  circumstellar 
material. 


The  Ko  program  is  based  on  non-local  thermodynamic  equilibrium  calculations  of 
ion  and  level  populations  employing  a  large  number  of  atomic  processes.  The  radiative 
transfer  for  both  line  and  continuum  radiation  is  treated  by  the  escape  probability 
formalism  (e.g.,  Krolick  and  McKee,  1978;  Drake  and  Ulrich,  1980).  This  method  works 
by  determining  self-consistent  optical  depths  and  level  populations  iteratively.  For  a 
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given  optical  depth  of  a  transition,  one  estimates  the  fraction  of  photons  that  escape  from 
the  gas.  Then  the  spontaneous  radiative  transition  probability  for  each  transition  is  set 
equal  to  the  Einstein  A  value  multiplied  by  the  escape  probability  and  the  rate  equations 
for  the  level  populations  are  solved.  If  the  optical  depths  from  the  newly  calculated  level 
populations  do  not  agree  with  those  originally  assumed,  then  the  calculation  is  redone 
with  the  new  optical  depth.  This  iterative  procedure  is  continued  until  self-consistency  is 
achieved. 

For  line  radiation,  there  are  three  probabilities  to  consider  in  computing  the  escape 
probability: 

1)   The  probability  that  the  photon  can  escape  without  taking  part  in  any  other 
atomic  processes,  J^/CO,  which  depends  on  the  optical  depth  of  the  line  and 
is  written,  using  the  notation  of  Ko  and  Kallman  (1994),  as 

WT/)  =  1  forx^lO"5,  (3-1) 

PesclJ(xl)  =  ^-^-  for  10~5  <x,  <l,  (3-2) 

2x 

and 

W*i )  = 7 '/i-T^t  f°r  T  " ]  (3-3) 

r  0.5Vlog(x,) 

T/VTI    1-2  +      — 

where  t„,  is  the  optical  depth  where  the  line  wings  become  optically  thick. 
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2)   The  probability  that  the  line  photon  is  destroyed  by  photo-absorption,  Pesclc , 

which  depends  on  the  line  opacity,  k,  ,  and  the  continuum  opacity  at  the  line 
energy,  kc  , ,  is  written  as 

/L,„=^^  (3-4) 


escl,c 


K/+Kc, 


where  kt  is  the  Thomson  scattering  opacity. 
3)   The  probability  that  the  photon  is  scattered  out  of  the  line  core  by  Thomson 
scattering  given  by 

escl.p  *  ^  ' 

The  escape  probability  for  a  line  photon  is  then  given  by 

For  the  continuum,  taking  the  line  analogy  for  the  escape  probability  of  the 

recombining  photons,  the  assumed  profile  function  of  the  recombining  photon  level  i  is 

q>,Jx)  =  cx-i  (3-8) 

where  jc  =  dzh  with  s  being  the  energy  of  the  photons  and  £,  the  ionizing  threshold  of  the 

level.  The  normalizing  constant  c  is  given  by 

$lJx)dx  =  \  (3-9) 

which  upon  integration  yields 

c  =  y-  1. 

At  the  absorption  edge  at  optical  depth  xa,  the  escape  probability  is 

\y,Jx)dx  (3-10) 

where  xc  is  defined  such  that 

which  yields 

*C  =  [(Y-1)VP  (3"12> 
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and  escape  probabilities 

^Cc(^)  =  [(Y-1K]('^  forxa>l/(Y-l)  (3-13) 

and 

^cCO  =  l  for  xa<  l/(y-l).  (3-14) 

For  recombining  continuum  the  profile  function  is 

T»M-     °»(X),X  (3-15) 

given  the  level  photoionization  cross  section  aioo  .  For  om(x)  equal  to  a  constant  times 
x"p, 

y=P  +  l  (3-16) 

which  yields 

^,/,(0  =  (PtJP/(M)  for   xa>l/P  (3-17) 

and 

C.(0  =  lfor  xfl<l/p  (3-18) 

with  P  =  3  for  hydrogenic  photoionization  cross  section. 

Heating  and  Cooling  Mechanisms 
The  Ko  program  was  designed  to  compute  the  emission  of  a  disk  in  an  X-ray 
binary.  In  these  systems  the  disk  surrounds  a  neutron  star,  and  the  accretion  of  disk 
material  onto  the  neutron  star  releases  X-radiation  which  is  then  absorbed  and  processed 
by  the  disk.  In  Algols,  the  gravitational  potential  well  of  the  primary  star  is  not  as  deep, 
and  X-radiation  is  not  produced.  Therefore,  certain  heating  and  cooling  mechanisms 
employed  by  the  Ko  program  are  negligible  when  applied  to  Algol  disks,  but  our  eventual 
goal  is  to  build  a  program  that  will  handle  a  wide  variety  of  matter- transferring  binaries 
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(including  X-ray  binaries  and  cataclysmic  variables)  and  the  Ko  program  is  well-suited 
for  our  work. 

The  viscous  interaction  in  the  circumstellar  gas  is  a  major  contributor  to  its  heating, 
and  the  method  of  computing  the  heating  was  described  in  Chapter  2.  The  Ko  program 
was  modified  to  read  the  viscous  heating  from  the  hydro  program,  add  it  to  the  other 
heating  rates,  and  compute  the  gas  temperature  for  thermal  equilibrium  (i.e.,  where  the 
heating  and  cooling  rates  are  equal).  Additional  heating  mechanisms  included  in  the  Ko 
program  are  photoionization  (due  to  radiation  from  the  primary  star)  from  all  levels  and 
Compton  heating,  which  is  significant  in  X-ray  binaries  but  not  in  Algols.  The  radiation 
field  of  the  primary  star  is  taken  to  be  that  of  a  blackbody  of  temperature  T  and  is  given 
by  the  Planck  function 


2s3      J_ 
nee       -1 


n*)=ih-^—  (3-i9) 


which  leads  to  the  photoionization  heating  rate 

element  ion  level  ^ 

*W=  IZI»/*L°^8)*  (3"20) 

k     J    , 
The  cooling  mechanisms  that  are  computed  are  radiative  recombination  cooling, 
line  cooling,  Compton  cooling,  and  bremsstrahlung  cooling,  the  latter  two  of  which  are 
very  small  in  the  Algol  regime.  The  radiative  recombination  cooling  rate  is 

element  ion  level  (   „     \  ._     (  0-c3 

*     j,        vhJ;,     Jct-  \  he  J 

and  the  line  cooling  rate  is 

clement   ion 

c,„,=  IIZI<%^.^  <3-22> 


k  m     j>i      i 
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where  e;/  is  the  line  energy  and  the  total  escape  probability  Pexcl  /(  is  given  by 

P'escl.^PescU+Pescl.pQ-PesCj)-  (3-23) 


CHAPTER  4 
RESULTS 


Spectroscopy 

Transient  Disk  Systems 

RW  Tau.  Kaitchuck  and  Honeycutt  (1982)  used  time  resolved  spectroscopy  during 
the  total  eclipse  of  the  short  period  Algol  RW  Tauri  and  found  that  the  disk  was  highly 
variable,  on  the  time  scale  of  one  orbital  period.  These  observations  dispelled  the  idea 
that  Algol-type  accretion  disks  are  axisymmetric  and  Keplerian  by  showing  that  the 
radius  of  the  disk  varied  from  1.1  to  1.7  times  the  radius  of  the  primary  and  that  line 
widths  were  at  least  a  factor  of  2  greater  than  would  be  expected  for  broadening  due  to 
the  rotation  of  a  Keplerian  disk.  Vesper  and  Honeycutt  (1993)  made  high  resolution  (~1 
A)  observations  of  the  Ha  line  over  the  entire  orbit  and  discussed  two  distinct  emission 
features  in  the  line  profiles.  One  emission  feature  ("Type  A")  is  strong  enough  to  be  seen 
in  raw  spectra  while  the  other  ("Type  B")  shows  up  only  when  the  light  of  the  two  stars  is 
subtracted  from  the  observed  spectrum.  The  subtraction  was  done  by  interpolating 
Kurucz  (1979)  atmospheres  for  the  primary  and  using  the  spectrum  of  HR  8857  for  the 
secondary,  and  then  convolving  them  with  limb  darkening  and  rotation  profiles. 

The  primary  of  RW  Tau  rotates  at  greater  than  twice  the  synchronous  rate  (Olson, 
1984;  Van  Hamme  and  Wilson,   1990).  Our  simulations  of  the  system  assume  the 

-8  -1 

parameters  given  in  Terrell,  et  al.  (1992),  and  a  mass  transfer  rate  of  10"  M0  yr  .  As 
Figure  4-1  shows,  the  stream  impacts  the  primary  star  and  rebounds  to  a  distance  of  about 
two  primary  radii  from  the  center  of  the  star.  As  mass  transfer  continues,  a  thin  ring  forms 
just  above  the  surface  of  the  primary  in  the  orbital  plane  (Figure  4-2). 
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Figure  4-1 .  The  stream  in  RW  Tauri  impacts  the  primary  star  at  phase  0.22. 
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Figure  4-2.  By  phase  0.32,  a  ring  of  material  has  formed  around  the  primary  star  as  matter 
continues  to  flow  from  the  secondary. 
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Figure  4-3  shows  the  Ha  line  in  RW  Tau  at  phase  0.222  as  observed  by  Vesper  and 

Honeycutt.  The  strongest  emission  feature  is  red-shifted  while  a  weaker  feature  is  blue- 
shifted,  and  spectra  taken  at  the  same  phase  during  other  epochs  show  similar  features. 
Near  phase  0.75  the  spectra  are  similar,  except  that  the  Doppler  shifts  are  reversed  (see 
Figure  3  of  Vesper  and  Honeycutt  (1993)).  At  phase  0.463  the  strongest  feature  is  blue- 
shifted  and  a  weaker  component  lies  near  the  rest  wavelength1  as  shown  in  Figure  4-4. 
Figure  4-5  shows  the  line  at  phase  0.685. 

For  the  computed  Ha  profiles,  we  chose  the  disk  as  shown  in  Figure  4-2  and  used 
the  Drake  program  to  compute  the  profiles  as  seen  from  different  orbital  phases.  It  must 
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Figure  4-3  The  Ha  difference  profile  for  RW  Tauri  at  phase  0.222.  The  dashed  vertical 
line  indicates  the  rest  wavelength  of  the  line.  (Data  courtesy  D.  Vesper) 


1  All  wavelengths  for  both  observed  and  computed  spectra  in  this  work  are  in  the  center 
of  mass  frame. 
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Figure  4-4  Same  as  Figure  4-3,  except  the  observation  was  made  at  phase  0.463. 
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Figure  4-5.  Same  as  Figure  4-3  with  the  observation  made  at  phase  0.685. 
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be  stated  that  we  have  not  attempted  to  formally  fit  the  observed  profiles.  At  this  point  we 

only  seek  to  show  that  the  computed  spectra  are  generally  similar  to  the  observed  ones. 

We  have  not  tried  to  fit  the  observations  by  varying  the  parameters  of  the  hydro  program. 

In  spite  of  this,  the  computed  profiles  do  match  the  observations  reasonably  well.  Figure 

4-6  shows  a  computed  spectrum  for  a  phase  0.463  view  of  the  disk  0.22  orbits  after  the 

onset  of  mass  transfer  (see  Figure  4-1).  In  the  observed  profile  we  see  three  peaks,  the 

strongest  of  which  is  blue-shifted  to  about  6554  A.  A  slightly  weaker  peak  lies  at  the  rest 

wavelength  and  a  much  weaker  peak  lies  at  6572  A.  The  computed  profile,  which  has 

been  smoothed  to  the  same  resolution  as  the  observed  one,  shows  similar  features  with  a 

peak  around  6552  A,  but  the  second  peak  is  slightly  red-shifted. 
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Figure  4-6.  The  computed  profile  for  RW  Tau  0.22  orbits  after  the  onset  of  mass  transfer 
and  as  seen  from  phase  0.463.  Figure  4-1  shows  the  disk  that  produces  this 
profile  and  Figure  4-4  shows  an  observed  profile  at  this  phase. 
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As  a  check  on  possibility  that  features  in  the  computed  spectrum  might  be  spurious 

and  dependent  on  the  number  of  gas  particles,  a  simulation  of  RW  Tau  was  run  with  the 

same  parameters,  except  that  the  number  of  particles  was  doubled.  Figure  4-7  shows  the 

computed  profile  for  the  resulting  mass  distribution  at  the  same  phase  as  the  profile 

shown  in  Figure  4-6.  Only  minor  differences  exist  between  the  two  profiles. 
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Figure  4-7.  Same  as  Figure  4-6  except  that  twice  and  many  gas  particles  were  used  in  the 
simulation. 


TX  Ursae  Majoris.  A  comprehensive  observational  study  of  the  circumstellar  gas  in 
this  system  was  done  by  Albright  and  Richards  (1993)  and  they  also  discuss  the  extensive 
observational  history  of  the  system.  Inspection  of  their  difference  profiles  shows  that  the 
circumstellar  emission  is  highly  variable,  with  rapid  changes  ocurring  over  intervals  as 
short  as  one  orbital  cycle.  Figure  4-8  shows  an  observed  Ha  profile  at  phase  0.306  and 
Figure  4-9  shows  the  profile  at  phase  0.631. 
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Figure  4-8.  An  observed  Ha  profile  of  TX  UMa  at  phase  0.306.  The  dashed  vertical  line 
indicates  the  rest  wavelength  (Data  courtesy  G.  Albright). 
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Figure  4-9.  An  observed  Ha  profile  of  TX  UMa  at  phase  0.63 1 . 
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Simulations  of  TX  UMa  were  done  assuming  a  mass  transfer  rate  of  10"   M0  yr" 

and  the  parameters  given  in  Terrell,  et  al.  (1992),  except  that  the  rotation  parameter  for 

the  primary  was  set  to  3.5  as  suggested  by  Albright  and  Richards  (1993).  Figure  4-10 

shows  the  distribution  of  circumstellar  matter  0.35  orbits  after  the  onset  of  mass  transfer. 


Figure  4-10.  View  of  TX  UMa  0.35  orbits  after  a  mass  transfer  event  has  begun. 

Ha  profiles  were  computed  for  the  disk  shown  in  Figure  4-10.  The  computed  profile  for 
phase  0.631  is  shown  in  Figure  4-11  and  has  some  features  in  common  with  observed 
profile  in  Figure  4-9  including  a  peak  at  6560  A  and  another  near  the  rest  wavelength. 
However,  the  strengths  of  the  two  peaks  are  reversed  and  the  computed  profile  shows  a 
small  peak  around  6570  A  that  is  not  seen  in  the  observed  profile. 
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Figure  4-11.  The  computed  Ha  profile  of  TX  UMa  at  phase  0.63 1  for  the  disk  shown  in 
Figure  4-9. 
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Figure  4-12.  The  computed  Ha  profile  at  phase  0.306  for  the  disk  in  TX  UMa  shown  in 
Figure  4-10. 
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U  Sagittae.  U  Sagittae  is  a  short-period  system  that  lies  in  the  transient  disk  region 

of  the  r-q  diagram  and  slightly  non-synchronous  rotation  with  F=1.31  (Van  Hamme  and 

Wilson,  1986).  The  Ha  difference  profiles  of  the  system  are  highly  variable  and  show 

strong  absorption  at  some  epochs  as  can  be  seen  in  Figures  4-13  through  4-17.  A  striking 

featue  of  the  observed  profiles  is  the  absorption  core  seen  at  some  phases. 

Simulations  of  the  system  were  run  with  the  parameters  in  Terrell,  et  al.  (1992). 

Unfortunately,  the  absorption  feature  makes  it  difficult  to  compare  the  computed  profile 

with  the  observed  ones.  Hoewever,  the  profile  computed  at  phase  0.499  does  share  the 

locations  of  peaks  and  overall  width  with  the  observed  profile  in  Figure  4-15. 
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Figure  4-13.  An  Ha  difference  profile  of  U  Sge  at  phase  0.171  during  a  time  of  high 
activity.  (All  U  Sge  data  courtesy  G.  Albright). 
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Figure  4-14.  Observed  Ha  profile  of  U  Sge  at  phase  0.219. 
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Figure  4-15.  Observed  Ha  profile  of  U  Sge  at  phase  0.499. 
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Figure  4-16.  Observed  Ha  profile  of  U  Sge  at  phase  0.672. 
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Figure  4-17.  Observed  Ha  profile  of  U  Sge  at  phase  0.888. 
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Figure  4-18.  The  computed  Ha  profile  of  U  Sge  at  phase  0.499  for  comparison  with  the 
observed  profile  in  Figure  4-15. 

Permanent  Disk  Systems 

AU   Monocerotis  is  a  long-period  (11.1    days)   system  that   lies  close  to  the 

permanent/transient  disk  dividing  line  of  the  r-q  diagram.  An  extensive  discussion  of  the 

system  is  given  by  Sahade  and  Ferrer  (1982).    Peters  and  Polidan  (1984)  analyzed  IUE 

spectra  of  the  system  and  found  evidence  for  a  high  temperature  accretion  region  with  an 

9  -3 

electron  temperature  of  approximately  1 00,000  K  and  electron  density  of  about  1 0  cm  . 
Figure  4-19  shows  an  Ha  profile  taken  by  G.  Peters  at  phase  0.22.  In  spectra  taken  about 
seven  months  earlier,  and  shown  in  Figures  4-20  and  4-21,  striking  absorption  features 
are  seen  at  the  core  of  the  profile 

Simulations  of  AU  Mon  were  done  using  the  system  parameters  given  in  Terrell,  et 
al.  (1992).  This  system  is  particularly  interesting  because  it  lies  very  close  to  the  Lubow- 
Shu  curve  in  the  r-q  diagram.  As  can  be  seen  in  Figure  4-23,  the  stream  barely  impacts 
the  star  and  forms  a  stable,  permanent  accretion  disk.  The  computed  profile  in  Figure  4- 
22  compares  well  with  the  observed  profile  of  Figure  4-18.  Of  course,  the  observed 
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Figure  4-19.  Ha  profile  of  AU  Mon  at  phase  0.22  with  the  flux  normalized  to  the 
continuum  (All  AU  Mon  data  courtesy  G.  Peters) . 
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Figure  4-20.  The  Ha  profile  of  AU  Mon  at  the  same  phase  (0.22)  as  the  profile  in  Figure 
4-18,  but  taken  about  seven  months  earlier. 


profiles  in  Figures  4-19  and  4-19  show  absorption  cores  and  cannot  be  modeled  by  our 
present  program  which  assumes  optically  thin  gas.  However,  the  width  of  the  computed 
Ha  profile  does  match  the  observed  ones  very  well. 
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Figure  4-21 .  The  Ha  profile  of  AU  Mon  at  phase  0.33,  taken  1  night  later  than  the  profile 
in  Figure  4-19. 
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Figure  4-22.  The  computed  Ha  profile  of  AU  Mon  at  phase  0.22  for  the  disk  formed  5.7 
orbits  after  mass  transfer  began. 
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Figure  4-23.  The  appearance  of  the  disk  of  AU  Mon  at  several  phases, 
a)  0.25  orbits;  b)  0.64  orbits;  c)  2.4  orbits;  d)  5.7  orbits. 

Polarization 


Interest  in  polarized  radiation  from  binary  stars  began  with  the  treatment  of 
Chandrasekhar  (1946),  who  showed  that  eclipses  would  break  the  disk  symmetry  and 
allow  limb  polarization  to  be  observed.  Other  workers  extended  the  Chandrasekhar 
treatment  (Code,  1950;  Collins,  1970,  1972)  and  predicted  very  small  eclipse  polarization 
effects.  Brown  et  al.  (1978)  treat  the  case  of  circumstellar  polarization,  and  derive 
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formulae  for  computing  the  Stokes  quantities.1  Wilson  and  Liou  (1993)  combined  the 
Chandrasekhar  and  Brown  et  al.  treatments  with  the  Wilson-Devinney  binary  star  light 
curve  model  (Wilson  and  Devinney,  1971;  Wilson,  1979,  Wilson,  1990),  and  were  able  to 
achieve  good  fits  to  the  Kemp  et  al.  (1983)  observations  of  Algol  by  assuming  fixed 
locations  of  circumstellar  material. 

The  hydro  program  has  been  interfaced  to  the  Wilson-Liou  polarization  program, 
making  it  possible  to  compute  polarization  curves  for  time-dependent  circumstellar  flows. 
Although  this  work  on  this  hydrodynamic  flow-polarization  (HFP)  program  is  at  a 
preliminary  stage,  the  results  demonstrate  that  polarization  observations  show  great 
promise  as  a  probe  of  circumstellar  flows.  The  ionization  treatment  at  present  is  rather 
simplistic.  The  ionization  is  computed  either  by  the  Saha  formula  or  by  assuming  that  a 
constant  fraction  of  the  gas  within  a  specified  radius  of  the  primary  star  is  ionized.  Future 
improvements  to  the  ionization  calculations  will  include  coupling  the  results  of  the  Ko 
program  for  electron  densities  to  the  polarization  program.  Figure  4-24  shows  the 
computed  polarization  curves  for  a  short  burst  event  in  SX  Cas,  assuming  Saha 
ionization.  A  longer  burst  event  results  in  a  dramatically  different  set  of  polarization 
curves,  as  seen  in  Figure  4-25.  An  even  longer  burst  event  is  shown  in  Figure  4-26.  The 
computed  curves  are  sensitive  to  the  assumed  ionization  scheme.  Figures  4-27  through  4- 
29  show  the  variation  of  the  computed  curves  for  the  constant  ionization  method  as  the 
radius  of  ionization  is  varied. 

Polarization  observations  of  Algol-type  binaries  are  relatively  rare  owing  to  the  fact 
that  necessary  combination  of  sensistive  polarimeters  and  large  aperture  telescopes  has 
not  been  developed.  Because  the  polarization  in  Algols  is  very  small,  large  numbers  of 


1  The  Stokes  quantities,  Q  and  U,  are  frequently  referred  to  as  the  Stokes  parameters, 
although  they  are  not  parameters.  We  reserve  the  term  "parameter"  for  quantities  of  the 
binary  that  are  estimated  from  observations  via  a  theoretical  model.  The  Stokes  Q  and  U 
are  directly  measured  quantities. 
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Figure  4-24.  The  polarization  curves  for  a  mass  transfer  event  in  SX  Cas  lasting  0.3 

orbital  periods.  The  upper  curve  is  the  U  curve  and  the  lower  one  is  the  Q 
curve. 


Phase 


Figure  4-25.  Polarization  curves  for  a  medium-length  (1.4  orbital  periods)  transfer  event 
in  SX  Cas. 
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Figure  4-26.  Polarization  curves  for  mass  transfer  event  in  SX  Cas  lasting  2.4  orbital 
periods. 
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Figure  4-27.  Polarization  curves  for  a  medium-length  (1.4  orbital  periods)  burst  event 
assuming  partial  ionization  of  all  circumstellar  gas. 
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Figure  4-28.  Same  as  Figure  4-4  except  that  the  ionization  occurs  within  a  radius  of  0.15 
of  the  primary. 


Figure  4-29.  Same  as  Figure  4-4  except  that  ionization  occurs  within  a  radius  of  0.1  of  the 
primary. 
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photons  must  be  collected  in  order  to  achieve  a  reasonable  signal-to-noise  ratio.  Although 

it  is  relatively  bright,  Algol  unfortunately  has  very  low  circumstellar  activity.  More  active 

systems  like  RW  Tau  are  too  dim  for  present  telescope/polarimeter  combinations. 

However,    interest   in   polarimtery   has   greatly   increased   in   recent   years,   and   the 

combination  of  sensitive  polarimeters  and  large  aperture  telescopes  like  the  10-meter 

Keck  Telescope  will  be  able  to  observe  polarization  in  many  of  the  more  active  Algols. 

Because  of  the  episodic  nature  of  these  gas  flows,  the  ideal  observing  platform  is  one  in 

space  where  weather  and  day-night  interruptions  are  eliminated  and  where  observations 

in  the  ultraviolet  can  be  made.  We  plan  to  further  develop  the  HFP  program  so  that  when 

the  observations  are  made,  we  will  have  a  tool  for  extracting  quantitative  information 

from  them. 

Conclusions 

The  goal  of  this  work  was  to  compute  the  spectral  radiation  arising  from 
circumstellar  gas  in  Algol-type  binaries.  Although  the  agreement  between  computed  and 
observed  spectra  is  not  perfect,  we  are  encouraged  by  the  similarities  seen  in  some  of 
them.  Of  course,  there  is  still  much  work  left  to  be  done,  and  the  work  presented  here  will 
serve  as  a  foundation  for  future  modelling  of  not  only  Algols,  but  other  types  of 
interacting  binaries  as  well.  The  first  step  in  using  a  new  model  is  to  become  familiar 
with  the  predicted  observables  for  a  wide  range  of  input  parameters  so  that  good 
estimates  of  the  input  parameters  can  be  made  before  beginning  an  impersonal  adjustment 
procedure. 

Some  parts  of  the  model,  in  particular,  the  radiative  transfer  program,  push  the 
limits  of  present  computers.  However,  further  coding  efficiencies,  combined  with 
inevitable  increases  in  computational  power,  will  make  it  possible  to  eliminate  some  of 
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the  simplifying  assumptions  necessary  in  this  work.  Enhancements  to  the  present  work 
will  be  made  for  many  years  to  come. 

In  this  project,  the  emission  lines  of  the  gas  were  compared  to  difference  profiles  of 
observed  stars,  where  the  light  of  the  two  stars  is  removed  from  the  observations.  This 
state  of  affairs  is  similar  to  the  way  light  curves  were  solved  with  the  Russell-Merrill 
model,  where  certain  effects  were  removed  from  the  observed  curve  before  solving  for 
the  system  parameters.  Unfortunately,  this  leads  us  in  the  case  of  our  emission  line  work 
to  view  the  stellar  light  as  a  contaminant  rather  than  a  source  of  information  on  the 
binary.  As  our  experience  with  physical  light  curve  models  has  taught  us,  it  is  better  to 
include  all  sources  of  light  in  the  model  and  predict  what  is  actually  observed  at  the 
telescope.  By  combining  our  work  with  the  Wilson-Devinney  binary  star  model  and 
Mukherjee's  (1993)  absorption  line  profile  model,  we  will  be  able  to  compute  observables 
that  can  be  directly  compared  to  observations.  This  combination  enables  us  to  better 
determine  the  parameters  of  the  binary  by  simultaneously  satisfying  different  kinds  of 
observations.  With  such  a  combined  model  we  could  simultaneously  fit  light  and  radial 
velocity  curves,  line  profiles,  polarization  curves,  and  X-ray  pulse  delay  curves.  The 
advantage  of  such  a  procedure  is  that  all  relevent  parameters  are  determined  in  a  logical 
and  self-consistent  way.  Another  advantage  is  that  simultaneous  solutions  may  extract 
information  not  readily  apparent  in  separate  solutions,  such  as  the  discovery  of  a  tidal  lag 
in  the  X-ray  binary  HD  77581  which  may  be  due  to  non-radial  pulsations  in  the 
supergiant  as  the  neutron  star  spirals  inward  toward  it  (Wilson  and  Terrell,  1994).  It  is 
hoped  that  the  present  work  will  spawn  new  efforts  both  theoretical  and  observational,  so 
that  a  coherent  picture  of  binary  star  evolution  may  be  achieved. 
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